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The problem ol soliton-soliton force is revisited. From the exact two solitons solution of a nonau- 
tonomous Gross-Pitaevskii equation, we derive a generalized formula for the mutual force between 
two sohtons. The force is given for arbitrary soUtons amplitude difference, relative speed, phase, 
■ and separation. The latter allows for the investigation of soliton molecules formation, dynamics, 

(N ■ and stability. We reveal the role of the time-dependent relative phase between the solitons in bind- 

^ , ing them in a soliton molecule. We derive its equilibrium bond length, spring constant, frequency, 

Q . effective mass, and binding energy of the molecule. We investigate the molecule's stability against 

' perturbations such as reflection from surfaces, scattering by barriers, and collisions with other soli- 

' tons. 
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I. INTRODUCTION 



' The old-new interest in the problem of soliton-soliton intertaction and soliton molecules has been increasingly 
^ ' accumulating particularly over the past few years. This is mainly motivated by the application of optical solitons as 
data carriers in optical fibers P, [2l and the realization of matter- wave solitons in Bose-Einstein condensates [4| . 
One major problem limiting the high-bit rate data transfer in optical fibers is the soliton-soliton interaction. On the 
"j^ . one hand, soliton-soliton interaction is considered as a problem since it may destroy information coded by solitons 
' sequences. On the other hand, it is part of the problem's solution, since the interaction between solitons leads to the 
I formation of stable soliton molecules which can be used as data carriers with larger "alphabet" @ . 
''O ' The interaction force between solitons was first studied by Karpman and Solov'ev using perturbation analysis 
Ch \ Gordon who used the exact two solitons solution [7|], and Anderson and Lisak who employed a variational approach 
O . Q. It was shown that the force of interaction decays exponentially with the separation between the solitons and 
I i ' depends on the phase difference between them such that in-phase solitons attract and out-of-phase solitons repel. 

This feature was demonstrated experimentally in matter- wave solitons of attractive Bose-Einstein condensates [j, 0] 
CO ■ where a variational approach accounted for this repulsion and showed that, in spite of the attractive interatomic 
^ , interaction, the phase difference between neighboring solitons indeed causes their repulsion Q. 

For shorter separations between the solitons, Malomed jlO] used a perturbation approach to show that stationary 
solutions in the form of bound states of two solitons are possible. However, detailed numerical analysis showed that 
. such bound states are unstable fllj . Stable bound states were then discovered by Akhmediev [l^ and a mechanism of 
creating robust three-dimensional soliton molecules was suggested by Crasovan et al. (iSl. Recently, soliton molecules 
were realized experimentally by Stratmann et al. in dispersion-managed optical fibers [5| and their phase structure 
' was also measured [l3|. Perurbative analysis was used to account theoretically for the binding mechanism and the 
molecule's main features [l^, [l^ . Quantization of the binding energy was also predicted numerically by Komarov et 
al. [13. InRefs.[ll[ll, a hamiltonian is constructed to describe the interaction dynamics of solitons. 

The mechanism by which the relative phase between the solitons leads to their force of interaction, and hence the 
binding mechanism, is understood only qualitatively as follows. For in-phase (out-of-phase) solitons, constructive 
?-H ' (destructive) interference takes place in the overlap region resulting in enhancement (reduction) in the intensity. As 
a result, the attractive intensity-dependent nonlinear interaction causes the solitons to attract (repel) [20| . A more 
quantitative description is given in Refs. [isl, [T6|. 

In view of its above-mentioned importance from the applications and fundamental physics point of views, we address 
here the problems of soliton-soliton interaction and soliton molecule formation using the exact two solitons solution. 
This approach has been long pioneered by Gordon [7] where he used the exact two solitons solution of the homogeneous 
nonlinear Schrodinger equation to derive a formula for the force of interaction between two solitons, namely 

A = -8exp(-A) cos(A(/)), (1) 

where A{t) is the solitons separation and A<j>{t) is their phase difference. This formula was derived in the limit of 
large solitons separation and for small difference in the center-of-mass speeds and intensities, which limits its validity 
to slow collisions. With appropriately constructed hamiltonian, Wu et al. have derived, essentially, a similar formula 
that gives the force between two identical solitons and reliefs the condition on slow collisions [l9| . 
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Here, we present a more comprehensive treatment where we derive the force between two sohtons for arbitrary soh- 
tons intensities, center-of-mass speeds, and separation. We also generahze Gordon's formula to inhomogeneous cases 
corresponding to matter-wave bright solitons in attractive Bose-Einstein condensates with time-dependent parabolic 
potentials [H, iS] and to optical solitons in graded- index waveguide amplifiers [2l|| . Many interesting situations can thus 
be investigated. This includes the various soliton-soliton collision regimes with arbitrary relative speeds, intensities, 
and phases. Most importantly, soliton-soliton interaction at short solitons separations will now be accounted for more 
quantitatively than before. Specifically, soliton molecule formation is clearly shown to arise from the time-dependence 
of the relative phase which plays the role of the restoring force. In this case, the force between the two solitons is 
shown to be composed of a part oscillating between attractive and repulsive, which arises from the relative phase, 
and an attractive part that arises from the nonlinear interaction. The time-dependence of the relative phase results 
in a natural oscillation of the molecule's bond length around an equilibrium value. The various features of the soliton 
molecule, including its equilibrium bond length, spring constant, frequency and amplitude of oscillation, and effective 
mass, will be derived in terms of the fundamental parameters of the solitons, namely their intensities and the nonlinear 
interaction strength. 

The two solitons solution is derived here using the Inverse Scattering method . Although the two solitons solution 
of the homogeneous nonlinear Schrodinger equation is readily known 0, US] ^ here we not only generalize this solution 
to inhomogeneous cases, but also present it in a new form that facilitates its analysis. The solution will be given in 
terms of the four fundamental parameters of each soliton, namely the initial amplitude, center-of-mass position and 
speed, and phase. The main features of the solution will be shown clearly such as the contribution of the nonlinear 
interaction to the actual separation and phase difference between solitons where it turns that the separation between 
the two solitons increases with logarithm of the difference between the amplitudes of the two solitons. Furthermore, 
the general statement that a state of two equal solitons with zero relative speed and finite separation does not exist as 
a stationary state for the homogeneous nonlinear Schrodinger equation, will be transparently and rigorously proved. 

Stability of soliton molecules is an important issue since, in real systems, perturbations caused by various sources 
such as losses, Raman scattering, higher order dispersion, and scattering from local impurities, tend to destroy the 
molecules. To investigate the stability of the soliton molecules described by our formalism, we have considered three 
situations. First, we studied the reflection of the molecule from a hard wall and a softer one. While for the hard wall 
the molecule preserves its molecular structure after reflection, it generally breaks up for the softer ones due to energy 
losses at the interface. Secondly, the scattering of the molecule by a potential barrier was also investigated. We show 
that the molecular structure is maintained only for some specific heights of the barrier. This suggests a quantization 
in the binding energy as predicted by Komarov et al. [17]. The oscillation period of the reflected molecule is noticed 
to be smaller than for the incident one. In addition, the outcome of scattering depends on the phase of the molecule's 
oscillation at the interface of the barrier. For instance, a dramatic change in the scattering outcome takes place if 
the coalescence point of the molecule lies exactly at the interface. In such a case, the otherwise totally transmitting 
molecule will now split into reflecting and tunneling solitons. Thirdly, we have considered the collision between a 
single soliton with a stationary soliton molecule. The effects of different initial speeds, amplitudes, and phases of the 
scatterer soliton were studied. It turns out that for slower collisions, it is easier for the scatterer soliton to break 
up the soliton molecule, while for fast collisions the scatterer soliton expels and then replaces one of the solitons in 
the molecule. The phase of the scatterer soliton plays also a crucial rule in preserving or breaking the bond of the 
molecule, which can be used as key tool to code or uncode data in the molecule. 

The rest of the paper is organized as follows. In section HTJ we use the Inverse- Scattering method to derive the 
two solitons solution of the inhomogeneous nonlinear Schrodinger equation and present the solution in the above- 
mentioned appealing form. The main features of the solution will be discussed in subsection III Al The center-of-mass 
positions and relative phases will be derived in subsections IIIBI and III CI respectively. The force between solitons 
will be derived in section III Dl where Gordon's formula will be extracted as a special case in subsection III D II and 
our more general formula will be derived in subsection HID 21 In subsection HID 31 we compare our formula with 
the numerical calculation. In section IIIII we show the possibility of forming soliton molecules, derive their main 
features in subsection IIII Al and investigate their stability in subsection IIII Bl We end in section ITVl with a summary 
of results and conclusions. The details of the derivation of the two solitons solution and the center-of-mass positions 
are relegated to Appendices A and B, respectively. 

II. THE EXACT TWO SOLITONS SOLUTION 

Matter-wave solitons of trapped Bose-Einstein condensates and optical solitons in optical fibers can be both de- 
scribed by the dimensionless Gross-Pitaevskii equation 

-I- + i {^{tf-^it)) xH{x,t)+g^e'<^'^\^{x,t)\H{x,t) = 0, (2) 
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where ^{t) is a dimensionless arbitrary real function. For matter-wave solitons, length is scaled to the characteristic 
length of the harmonic potential, ax = i/S/mWa,, time to 1/uJx, and the wavefunction to Xj ^J^a^j^ujijuix^ 

where uj^ and a;_L are the characteristic frequencies of the quasi one-dimensional (w^ 3> Wj,) trapping potential in the 
axial and radial directions, respectively. In these units, the strength of the interatomic interaction will be given by 
the ratio go = Os/a^, where is the s-wave scattering length. For the case of optical solitons, the function \['(a;,i) 
represents the beam envelope, t is the propagation distance, x is the radial direction, and the intensity-dependent 
term represents the Kerr nonlinearity. In this case, scaling is in terms of the characteristic parameters of the fiber, 
as for instance in Ref. [19]. The specific form of the prefactors of the inhomogeneous and nonlinear terms guarantees 
the integrability of this equation |23l [23 |. For the special case of 7(t) — 0, the homogeneous case is retrieved. Other 
interesting special cases have also been considered |23l425| . 

As outlined in Appendix \^ we use the Darboux transformation method to derive the two solitons solution of this 
Gross-Pitaevskii equation, which can be put in the form 

^{x^t) = y sech[an - x,^^{t))\ 



+ 



n2a22(O ^i(0oi+0O2+02(3:,t)+tan-^(ga)+tan-i(gA)) 

T C i J 
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a22{t){x - Xcrn2{t)) + - log 



2 \ 1 
2 



-I- a/ 



(3) 



where 



J = 1,2, 

'^cmj (0 — -^cmj(07 

<l>j{x,t) = y{t) (4e-270^2 ^2^2) _ 1^^_^^.(^)2^(^) _^ V,„,i{t){x - X,„,j{t)) - ^{x - X,^i{t))mt), 

a\ = ,fi + e^" cos z, 
"2 = /2 + e^" sinz, 

"4= 2^ - '^^idoe^" sinz, 

fi = ("-2 + ni)go + (712 - 7ii)goe^, 

/2 = 2(r;2-?7i)e"^" (1 + e^), 

/a = -(?^2 - ?^i).go - ("2 + ni)goey, 



Vn 



e'^'^5o(("i - "2)2: - («ia;cmi(i) ~ "2a:cm2(i))), 



go(("i + n2)x - (nia:cmi(i) + "-2a:cm2(i))), 



Z = -002 + Zii - Z22, 

Vj = Vj +Xji{0), 
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The solution is put in this suggestive form to facilitate its analysis. The first sech part corresponds to the 
exact single soliton solution with center-of-mass position Xcmi{t), width l/aii(t), phase ^oi + 4'i{^^t)i S'Hd nor- 
malization ni. Hence, xi and vi correspond to the initial center-of-mass position and speed, respectively. The 
second sech term contains the same features in addition to a shift in both the center-of-mass position and phase. It 
should be noted, however, that Xcmjit), Vcmjit), rij, and (j)j{x,t) correspond to the ccntcr-of-mass position and speed, 
normalization, and phase of the single noninteracting solitons. Due to the interaction between solitons, these four 
characteristic quantities may not correspond exactly to the values of the same physical quantities as they did for the 
single soliton solution. For instance, Xcmiit) will not correspond to the center-of-mass of one of the solitons. Instead, 
the soliton may be shifted from that position due to the interaction with the other soliton. In the following, we 
present a detailed analysis of the locations and phases of the two solitons. 

A. Main features of the solution 

Inspection shows that there are two main regimes for the two solitons solution, namely the regime of resolved 
solitons and the regime of overlapping solitons. In the former case, the center-of-mass concept is well defined and 
analysis of the relative dynamics becomes feasible. The solitons are considered resolved as long as the two main peaks 
are not overlapping, which means that partial overlap may occur in this regime. The analysis in this section assumes 
the resolved solitons regime. 

In the resolved solitons regime, the argument of the second sech term of Eq. ([3]), namely q = 022 (^) (2: — 2:cm2(^)) + 
(1/2) log [{a1 + Q;2)/(ck§ -I- , simplifies to a function with three roots. The fact that the sech function is peaked 
at the roots of its argument, leads to that the second sech term corresponds to three "solitons". This is shown in 
Fig[TJ We denote these solitons as the "left", "central", and "right" solitons with peak locations at xi, Xc~ Xcrai{t), 
and Xr: respectively. We notice that the central soliton is located at the position of the soliton of the first sech term, 
namely near Xc-nxi{i)- Further inspection shows that the two solitons at this location are out of phase and interfere 
destructively such that they do not appear in the total profile. Therefore, the two solitons that our solution of Eq. ([3]) 
describes are in fact the left and right solitons arising from the second sech term. This is different from what one 
aught to conclude from the form of the exact solution, namely that the first sech term corresponds to one soliton and 
the second sech term corresponds to the other soliton. 

In general, the center-of-mass locations of the left and right solitons, xi and Xr, do not match Xcmi(^) and Xcm2{t)- 
Typically, xi will be shifted to the left of Xcmi(i), Xr will be shifted to the right of Xcm2(i), while Xc remains near 
Xcmi{t)- The amount of shift depends mainly on the normalization difference ^2 — "-i and the relative speed V2 — vi, 
as will be shown in the next subsection. An interesting general and exact result is that, for the homogeneous case 
7(t) = 0, a state of two equal solitons, ni =11,2, with zero relative speed, vi — V2, and finite separation, does not exist 
as an exact solution of the Gross-pitaevskii equation. This can be proven by substituting ni = 712 and vi = V2 in q 
to find that the right and left solitons migrate to 00 and —00, respectively, while the center-of-mass of the central 
soliton matches exactly Xcmiit). Furthermore, we show in section fll CI that the phase difference between this central 
soliton and that of the first sech term equals, in this case, tt; guaranteeing their destructive interference. Thus, the 
three solitons disappear in such a special case and \E'(a;,t) becomes the trivial solution. 

In view of the above, a need arises to derive formulae for the center-of-mass positions, xi and Xr in terms of the 
solitons parameters, which will then be used to derive the force between the two solitons. 

B. Center-of-mass positions 

In this section, we derive formulae for the three roots of q which correspond to the locations of the left, cen- 
tral, and right solitons, xi, Xc, and Xr, respectively. To facilitate the derivation, we define X — exp (y™) — 
exp [ie''(*).go {{ni - ^2)2; - (ni Xcmi{t) ~ »^2 Xcm2{t)))] , Y = exp (y) = exp [e^(*)ni5'o(a; - Xcmiit))] , Ud = n2 - ni, 
Us = ni + n2, rjd = 'r]2 — ^i, and r/s = J?! + '72- The equation q{X^ y) = is a third-order polynomial in both X and 
Y . In principle, this equation can be solved algebraically for X or Y . However, extracting x from the resulting three 
roots will not be possible analytically for general ni ^ n2- Alternatively, and as can be seen from Fig. [TJ we can 
exploit the simple linear behavior of q near its roots. 

Assuming, without loss of generality, Xcm2{t) > a;cmi(^): ^^d noting that in the resolved solitons regime the solitons 
separation \xcm2{t) — Xcmi{t)\ is large, we argue in Appendix IB II that X ^ 1 for all x, F ^ 1 for a; ~ Xr, F ~ 1 for 
X ^ Xc, and Y <^ 1 for a; ~ x;. Based on this, the center-of-mass position Xr can be derived from a Taylor expansion 
of q in powers of large X and Y . Expanding q in powers of large X only and leaving Y arbitrary, accounts for Xc and 
xi simultaneously. Keeping terms up to first order of 1/X, we show in Appendix IB II that the position of the right 
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soliton is given by 



^ =^ Jt) 2e-^'*)logy. ^ 4y.":+"- (gon.e^°cos^.-2r?,sinz.)e-^W+-^° ^_i,„,^(„^_,^^,.w 
go{nd + n,) go{nd + ris) {g^nd^e^-^o + Atj^^) 



(4) 



where 



and the position of the the left sohton is given by 



where = a;cm2(^) — a;cmi(i), = z(xr)^ and C+ and are given in Appendix IB II 

The second and third terms on the right hand side of Eqs. ^ and ^ account for the shift in the center-of-mass 
position with respect to the single soliton ones. The third terms are much smaller than the second ones since they 
decay exponentially with the solitons distance Xd- In the limit — > and r]d — > 0, both logF^ and logy+ take the 
form log [(n^ + 4e~2'^''?7^/(7Q)/n^]. Thus, it is obvious that for = ?yd = 0, = oo and xi = — oo. This agrees 
with our earlier result that two equal solitons with zero relative speed and finite separation, do not exist as an exact 
solution. 



C. Relative Phases 



In this section, we calculate the phases of the left, central, and right solitons in reference to the phase of the soliton 
of the first sech part in the two solitons solution. For simplicity, the special case of ^ 1 and rjd 1 will be 
assumed. The phase difference between the left, central, and right solitons on one hand and the soliton of the first 
sech term on the other hand is generally given by 

A0 = (/.2(a;,i) +<?^02 +tan-i(— ) +tan-i( — ) - ipiix.t), (7) 

ai as 

which by observing that for all x 

(f>2{x,t) + (j)02 ~ (f>iix,t) ^ z (8) 

reduces to 

A0 = z + tan-i( — )+tan-i( — ). (9) 
ai as 

This expression gives the phases of the right soliton (f>r, the central soliton 0c, and the left soliton (f>i, for 
and xi, respectively. 

To calculate these phases we express the parameters Q!i_4 in terms of X and Y, as follows 

Oil —risgo + ndgoY + X cos z, (10) 



a2=Vde-""'{l + Y) + X sinz, 



Ud + UsY Y 
as = (ris - rid) go -rp cos z, 



Us - Ud 



X 



(11) 
(12) 



2r,de-'yo{l + Y) Y . 

aA^—, r {fis - Ud) go— s\nz. (13) 

(lis -rid) go X 

In general, X > 1 for ah a;, but F > 1 only for x > Xcmi(t), and Y = X e"'''*'' 9o^^'^^-^"^2{t))/2 > x for x > Xcm2{t), 
as shown in Appendix IB II 
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I. Phase of the right soliton (/i^: 

In this case x > a:cm2(i) which leads to Y ^ X and thus ai = UdgoY, «2 = '^Vd Y , = —nsY/iris — rid), 
a4 = 2'qdY e~^° / {iris - nd)go)- Therefore 



tan-i ^ = tan 

Ql 



^1 / 2-nie—<0 \ 
V go rid J 



0, 
7r/2, 
0[nd,r]d), 



nd 7^ 0, rid = 0, 
nd ^0,r]dy^ 0, 
nd y^O,r]dj^ 0, 



and tan"^ ^ = tan"^ I ^jjd^Jl] = which gives 

0r==7r + 6l-z, (14) 
where is a function that depends on the ratio rjd/nd for nonzero rjd and rid- 
II. Phase of the central soUton 4>c- 

In this case x — Xcmi{t) which gives Y = 1, and hence tan^^ ^ = z, and tan^^ ^ — tan^^(— 1,0) — tt. 
Finally, we get 

0c = TT. (15) 

The last result shows that, for rjd = nd — 0, i.e., two equal solitons with zero relative speed, the central soliton and 
the soliton of the first sech term of the two solitons solution are out of phase and therefore interfere destructively. 



III. Phase of the left soliton 

In this case x < Xcmi{t) which results in K ^ 1 and ai = X cosz, a2 = X sinz, — ^nd/{ng — nd), 
a4 2rid e-^o/{ns - nd) go- Therefore, tan"! |^ = z and tan"! ^ = tan'^ ( 2r;rfe~-"o/go ) which gives 

(j)i^2e-z. (16) 

The phase difference between the left and right solitons is thus given by Acf) — (t)r — (t>i = 29 — z. Noting that 9 — 
for rjd — and small but finite rid, and 9 = tt/2 for rid — a-nd small but finite rjd, we finally conclude that the phase 
difference between the two solitons is given by 

{TT - z, nd = 0,rid^ 0, 

(17) 
-z, r]d = 0,nd7^ 0. 

The above results are verified in Fig. [2l where we plot A(^, (jjr, 4>i, and (pc versus x. The agreement between our 
estimated values and the exact curve is evident. 
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D. Soliton-soliton force 

In this section, we use the results of the previous two subsections to derive the force between the left and right 
solitons. The force is proportional to the acceleration of the solitons separation 



A = — Xi 

= Xd(t) 



ae 



-lit) 



+ COs{Zr{t))e-^3onaxait)e-'^''>-,(t)+2-to 

+ P2 cos{zi{t))e-isonaxait)e-'^*^~'rit) 
+ 133 sin(z,(t))e"^S''"''^''(*)'=""*'-'^(*)-'^° 



where 



2{{n, - Ud) log(y.) + {rid + n.) log(y+)) 

and the coefficients /3i_4 are given in Appendix IB 21 The first two terms on the right hand side of Eq. (|18p are the 
dominant ones since they correspond to the noninteracting solitons separation, Xd (t) , and their logarithmic shifts oc a 
arising from the interaction between solitons. The time-dependence of A originates from j(t), Zr^i{t), and Xd{t). The 
acceleration can thus be derived 

m = iiitf-lit)) m + ^ [/35sin(zKi)) +/36sin(z,(i)) +/37COs(z,W) +/38COs(z.(i))] 

(20) 

where the coefficients (3^-$ are given in Appendix IB 21 In the last equation, we have used Eq. with only the 
first two terms of its right hand side to substitute for Xd(t) in terms of A(t) in the exponential factor. The first 
term on the right hand side of Eq. (l20l) corresponds to the force due to the external potential which vanishes for the 
homogeneous case. The rest of the terms correspond to the force of interaction between the solitons. The interaction 
force depends, as expected, on the phase difference of the two solitons and decays exponentially with their separation. 
It should be noted that this equation is a generalization of Gordon's formula in two aspects. First, it is derived 
for a time-dependent inhomogeneous medium. Secondly we have, essentially, no restriction on the difference between 
the two solitons amplitudes and speeds; apart from some extreme cases which were mentioned in Appendix IB II and 
will be discussed further below. 



1. Gordon's formula 

For the homogeneous case, 7(i) = 0, and in the limits Ud — > and ijd — > 0, the acceleration formula, Eq. (|20p . 
simplifies considerably. An apparent inconsistency occurs when switching the order of these two limits, namely 

lim lim A(<) = (n^5o)^e^3"»9°'^ cosz, (21) 

while 

lim lim A(t) = i (n,5o)^e-3"»9''^ cosz, (22) 

which differ by an overall minus sign. The conflict is resolved by invoking Eq. (|17p where it is shown that in the first 
case: z = A0, while in the second case: z = tt — Acf). Therefore, the two approaches agree on the following result 

Ait)^-^{n,gofe-i"^<'°^ cosAq^. (23) 

This is essentially Gordon's result, Eq. ([T]), since in his derivation Gordon took go = 1 and soliton amplitude rij^fg^jl — 
1, j = 1,2, as can be seen in Eqs.(l) and (6) of Ref. 0. Substituting go = 1 and = 4 in the last equation, it 
becomes identical to Eq. ([T]) . It should be mentioned here that Eq. ([IJ was also derived in Ref. 6] using a perturbation 
analysis based on the Inverse Scattering method, and in Ref. Q using a variational calculation. 
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2. Our formula 

For nonzero ^{t) and in the limits — )• and rjd — > 0, the acceleration formula takes the form 

K{t) = (7(t)2 - A{t) - ^ (go n,)^ g-hon^A^t) e-'W+37(t)-70 + g7o ) cos(A0). (24) 

This is a generalization to Gordon's formula for the inhomogeneous case as modeled by Eq. ([2]). Depending on the 
specific form of 7(i), the two force terms, namely the external (first term) and the interaction (second term), can be 
repulsive, attractive, or oscillatory. In addition, the phase difference A0(i) also depends on j{t). It is established 
in the homogeneous case, as will also be shown in section IIIIl that the time-dependence of the phase difference 
is responsible for binding the two solitons in the soliton molecule. Here, in addition to the possibility of forming 
soliton molecules, the dependence of A(/)(t) on 7(i) allows for controlling the parameters of the molecule such as its 
equilibrium bond length, period, and spring constant. This and possibly other interesting phenomena will be left for 
future investigation. 

3. Comparison with numerical calculations 

To obtain an estimate of the accuracy of the general acceleration formula, Eq. (|20|) . we calculate numerically the 
acceleration from the exact two solitons solution, Eq. ([3]), and compare the two results. The distance between the 
two solitons of the function |5'(x,t)p is determined using a numerical algorithm that employs our formulae for xi and 
Xr given by Eqs. and © to calculate seed values. The distance is then differentiated numerically twice at t — 0. 
In Fig. [3l we compare the two results. Good agreement is obtained for \rid\ < 1. The analytical solution diverges at 
rjd ~ ±1. The value at which divergence takes place is set by the specific choice of parameters in Fig. |31 As pointed 
out in section fll Bl the divergence occurs due to the merging of the central and left solitons. This artifact divergency 
can be remedied by associating the location of the local maximum of q to xi once this maximum has reached the 
a;-axis from above. 

Restricting our study to the region where agreement is obtained, we interestingly notice that the acceleration is 
oscillating between positive and negative values. This means that the force between the solitons is oscillating between 
attractive and repulsive. The possibility of having attractive forces for finite rjd is particularly interesting; for two 
solitons with nonzero relative positive speed, i.e. the solitons are initially diverging from each other, the force between 
them is attractive. This suggests that, if the force remains attractive for sufficient time, the two solitons will slow 
down and eventually converge at some point. If true, this should occur at small distances since the force decays 
exponentially with distance and when the two solitons are allowed to diverge even for a short while, the force might 
be weakened such that the two solitons can not return back. To be able to judge on such a possibility, we need to 
know what happens to the acceleration at later times. To that end, we calculate numerically the acceleration in terms 
of rjd and t. The result is plotted in Fig. |31 where it is clear that the acceleration indeed decays with time for all rjd- 
This leads to that any nontrivial effect of the oscillating force is most likely to take place at short solitons separations. 
This is what we find in the next section where the possibility of forming stable soliton molecules is pointed out. 

III. SOLITON MOLECULE 

We have shown in the previous section that, as a result of the solitons time-dependent relative phase, the force of 
interaction between solitons is oscillating between repulsive and attractive. Since the force decays exponentially with 
the solitons separations, this oscillation will have a tangible effect only when the two solitons are close to each other. 
In this section, we investigate the force of interaction between solitons for short solitons separation. In such a special 
case, Eq. (1181) takes a simple form that accounts for the solitons separation in terms of their relative phase. Using 
this formula, we show that the solitons will be bound to oscillate around some equilibrium distance where the phase 
plays the role of the restoring force. Comparison with exact numerical calculations shows that this formula is accurate 
for almost the full range of the solitons separation, except at the coalescence point (if any) . In subsection IIII A[ we 
discuss the main features of the resulting soliton molecules, and in subsection IIII Bl we investigate numerically their 
stability in different scattering regimes. 

To focus on the role of relative phase, we simplify the analysis by restricting our treatment to the homogeneous 
case, j{t) = 0, and zero relative speed, rid = Vd = V2 — vi = 0. We also set xi = xi so that any separation between 
the solitons to be as small as possible which in this case arises only from the logarithmic shifts (~ a in Eq. (US])). In 
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this case, the sohtons separation A, given by Eq. (|18p . simphfies in the hmit ^ to 



gang 



■ los 



1 + 2 go cos - 30 nrf ris n + (n^ 50) 



4 .log^^°"' 



and the acceleration is given by 



A(t) 



1 ^ ^3 2 Tig go + {us go)^ cos (g g§ rig t) + cos (| gg t) 
8 1 + 2?is go cos (ig^nrfn^i) + (nsgo)2 



(25) 



(26) 



It is noted that for go 3> 1 or go ^ 1, Gordon's formula is retrieved, but here with an explicit time-dependence 
of the phase. At/) = Ud Ug go t/^- This acceleration formula deviates considerably from Gordon's formula for go ~ 1. 
Specifically, for go = 1, A diverges to —00 at cos (n^ go t/8) = — 1, which indicates that the two solitons coalesce. 
This is confirmed below by examining the exact solution at this condition. We note here that an approximate 
expression for the solitons separation was also derived in Refs. d, [l^]. In addition, our predicted molecule's 
oscillation frequency (see Eq. (p3)) below) agrees with these references. 

To verify this feature, we calculate numerically the distance between the two solitons directly from the exact 
solution, Eq. ([3]), for different values of n^go. For ri^ go — 2.5, the density plot in Fig. [S^, shows a soliton molecule 
of two clearly resolved solitons with a separation oscillating around some nonzero equilibrium distance. Approaching 
the solitons coalescence point with n^go — l-Sj the density plot in Fig. |6^, shows the two solitons approaching each 
other more than the previous case. Furthermore, this figure shows a slight bounce back by one of the solitons in the 
region of collision. Approaching further the coalescence condition with ri^ go — 1.25, we indeed observe in Fig. [7^ 
that the two solitons merge almost completely. For a more quantitative comparison, we calculate numerically the 
center-of-mass trajectories of the two solitons. We show the trajectory curves in the density plots of Figs. [5^-[7^- 
In Figs. [S}d-[7)d, we plot the solitons separation obtained from formula (1^51) and the numerical trajectories obtained 
from the exact solution. It is clear from these figures that this formula agrees well with the exact soliton separation 
except near the collision region. In Fig. [5)d, the two solitons remain away from each other during the collision, and 
therefore good agreement is obtained with the exact result even in the collision region. In Fig. |6}d the two solitons 
approach each other further such that formula does not account for the above-mentioned slight bounce of one 
of the solitons. In Fig. [71d, agreement with the exact solution in the collision region is qualitative. We found that at 
the condition ris go = 1 and for rid <C n^, the analytic curve overlaps with that of the exact solution; apart from the 
horizontal segments where formula (I25p diverges to —00. Further insight is obtained by plotting the density profile of 
the soliton molecule at some specific times, as shown in Figs. [S};-!?!:. In Fig. [St, we observe that the initial amplitude 
imbalance is never removed during the dynamics. Instead, it becomes maximum when the two solitons are closest to 
each other. In addition, we notice that the oscillation amplitude of the larger soliton around its equilibrium position 
is larger. Figure [Ht shows clearly the soliton bounce which takes place in the time interval t = hh to t = 80. In 
these subfigures we plot two vertical dashed lines that indicate the position of the solitons at the closest approach. It 
is clear that after first closest approach at t = 55, the right soliton bounces back with a maximum displacement at 
t = 66. In Fig. [T];, it is shown that, although the two solitons coalesce, two small symmetric wings appear. A detailed 
examination of these wings shows that they are the remnants of the two solitons after they coalesce and they both 
bounce back in the collision region similar to the case of Fig. [51 

It is also instructive to show the dynamics of the phase profile during the molecule's oscillation. This is shown 
with the contour plots of Fig. |8] which correspond to the molecules of Figs. [5][7l In Fig. [8^, which corresponds to 
Fig. m the two solitons start initially in phase. By time the phase of the right soliton, which is the one with higher 
intensity amplitude and larger oscillation displacement, starts to exceed that of the left soliton. At the point of closest 
approach, the phase difference is exactly tt. After that point, the two solitons diverge again, the phase difference 
starts to decrease, and the cycle is repeated. Similar behavior is seen in Fig. [8]d. However, in Fig. [S}:, where the two 
solitons coalesce for a considerable amount of time, the phase difference during the coalescence time is zero. It is thus 
not completely understood why, in this case, the two solitons still repel each other and eventually split. 



A. Molecule formation and dynamics 



Having established the existence of the soliton molecule from the exact two solitons solution and derived a formula 
that describes its bond length, here we use this formula to examine more closely the propertie of the soliton molecule 
and its mechanism of binding. 

It is clear from Eq. (pS)) that the sinusoidal time-dependence of the solitons relative phase leads to a force of 
interaction that oscillates between attractive and repulsive and hence allowing for soliton molecule formation. Further 
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details of the mechanism of binding will be uncovered by expressing the acceleration, A(t) in terms of A{t) by 
substituting for cos {ug n^t^Q t/8) from Eq. (^51) into Eq. (^5]) . to get 

A = l.gon. (1 - glnl) ' (^^) ' e-^on^A _ _^ 2 (^2^2 ^ ^-i,„„.A^ (27) 

This shows that the interaction force between the two solitons is the resultant of an attractive part and a repulsive 
part. The equilibrium bond length, defined by A(A = Aeq) = 0, is given by 



= log 2 t 2 2 ■ 28 



In consistence with our previous result, the equilibrium bond length diverges as — logn^. Solving the last equation 
for and then substituting in Eq. ([27| . A simplifies to 



-39on=(A-%l)^ 2 _ p-igon.A 

For small amplitude oscillations, A ~ Acq, the last equation gives 



(29) 



A = -^alnl {glnl + l) (A - Aeq)e-3^/on.A.,^ (30) 

The restoring force (~ A) originates from the phase-dependent terms, cos (n^ rig */8). This appealing form of the 
force of interaction shows that the force between the solitons is of Hooke's law type with a spring constant 

^=S5o'-^.9o"' + l)e"^^""=''"^ (31) 
where m is the bare mass of the molecule. Expressed in terms of Ug and n^, the spring constant takes the form 

, _gW,nl{g',nUl)' 



64{g^,nj-lY 



■m, (32) 



which shows that fc = for rid = 0; corresponding to a soliton molecule of infinite bond length. Furthermore, k 
diverges for ng go — 1, which signifies soliton coalescence, as we have pointed out in the previous subsection. Since 
the frequency of the soliton molecule is given by 

uj = ^ndnggl, (33) 
and the spring constant is given in Eq. p2p . the effective mass m* = k/ui'^ will be given by 

(34) 
(3o'^?-ir ' 

which again diverges at the soliton coalescence condition, Ug go = 1. Having determined the main properties of the 
soliton molecule, we can now return back to Eq. to express A as 

A = A„ + ^log f 9'o-'s+2gongcosi.t) + l\ ^^^^ 
gong \ [gorig + 1) ^ / 

where Aq is the initial value of A, which is given by the solitons parameters through 

Ao = A.q + ^logf-l±^y (36) 
gong \{l-gongyj 

The amplitude of the oscillation A^ax — Aq — Acq is thus given by 

A.. = ^ log f -^^^l , (37) 
gong \{l-gongyl 
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which gives an elastic potential energy 



E=7^k A„,^^ = - go -; :t-^ log ^2 ■ (38) 



r, "max a JO '"d 1 i 2„ 2 / & /i \2 

^ f> \i -50*^5 / y(l - gaUs) 

It should be noted here that this is equal to the mechanical energy since the initial speed vanishes, A(0) — 0. The 
fact that the potential energy diverges at the coalescence condition 50 = 1 is a gain an artifact of the calculation, 
but it at least indicates that the bond is tighter than cases where rtg go 3> 1 or Us go ^ 1. 



B. Stability 

Here, we investigate the stability of the soliton molecule against break up in the following three collision regimes: 
i) reflection by a hard wall, ii) crossing a finite potential barrier, and iii) collision with a single soliton. To that end 
we solve the Gross-Pitaevskii equation, Eq. ([2|), numerically. As an initial state, we use, for cases i) and ii) the two 
solitons solution, Eq. ([3]), which represents the soliton molecule. For case iii), we use the superposition of the exact 
single soliton, Eq. (jA4l) . with the two solitons solution. 

Before starting the discussion of results, we point out that in Figs. [MTU we present the results of this section 
using spacio-temporal density plots. Since the solitons are too thin compared to the spacial range that we consider, 
a density plot with full spacial and time ranges will not show a clear solitons peak density or center-of-mass path, as 
Fig.O; shows. To solve this problem, we restricted the density plotting to a finite range of |5'(a;,i)p, namely between 
0.025 and 0.15 corresponding to the upper part of the solitons peaks. This results in an easier tracking of both the 
solitons peak density and center-of-mass path, as shown in Fig. [9^,b,d and the rest of subsequent figures. 

For reflection from a hard wall we solve Eq. © with a potential step of the form 

{Vb, x<xq, 
(39) 
0, x>xo, 

where Vq and xq are the hight and location of the potential wall, respectively. The result of reflection from this hard 
wall, with Vo = 100, is shown in Fig. ^1^. The soliton molecule preserves its molecular structure but with different 
characteristics. The solitons in the reflected molecule do not coalesce as in the incident molecule. In other words, 
the equilibrium bond length becomes larger. The density plot shows that initially the two solitons are of comparable 
intensities. After reflection, the brighter color of the left soliton and darker color of the right soliton indicate that the 
left soliton acquires higher intensity on the expense of the right soliton. We also notice that the left soliton performs 
two reflections from the potential interface. After the first reflection, it collides with the right soliton and then collides 
with the potential interface for the second time. 

The picture becomes different when the hight of the wall is reduced to Vb = 0.075, as shown in Fig.[5]3. The soliton 
molecule breaks up after reflection. This is due to loss of energy at the interface of the potential. Part of the soliton 
molecule transmits as a nonsolitonic pulse that broadens and decays in intensity by time. By plotting |\E'(2:,i)p in 
Fig.Ht with its full range, we can see the nonsolitonic part as the left- and right-going two red ejections corresponding 
to the transmitted and reflected nonsolitonic pulses, respectively. In Fig. [SJi, we combine Fig. ^jp and Fig. [SJ: to show 
the locations of the nonsolitonic ejections with respect to the solitons centers. In the case of reflection from a hard 
wall, the nonsolitonic ejections are essentially not present which results in the stability of the molecular structure. 

For reflection from a potential barrier we solve Eq. ^ with the potential 

{Vb, Xq - d < X < Xo, 

(40) 
elsewhere, 

where Vq, d, and xq, are the height, width, and location of the right side of the barrier, respectively. In Fig. [TUl we 
show the many different possibilities that result when the height of the barrier is changed. The free evolution case 
with Vb = is shown as a reference plot. The full reflection case is shown for Vb = 100, which is similar to the previous 
case of reflection from a hard wall. Reducing the height of the barrier to Vb = li we notice that the soliton molecule 
breaks up after reflection. As pointed above, this is due to the nonsolitonic ejections taking place at the interfaces 
of the potential. Reducing the height of the barrier to Vb = 0.5, a sign of solitons recombining appears in the form 
of a soliton molecule of a short lifetime. At Vb = 0.465 a stable molecule is remarkably formed with a considerably 
shorter period than for the incident molecule. We have confirmed numerically that this molecule remains stable for 
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much longer time provided that the soliton molecule remains sufficiently far from the boundaries of the spacial grid. 
This unique structure remains for some small domain around Vq = 0.465, but is lost for Vq — 0.45, where the molecule 
breaks for long evolution times. Decreasing the height of the barrier to Vq = 0.4, the molecule breaks at the interface 
and splits into a reflected and transmitted solitons. For Vq = 0.075, the molecule breaks at the interface, but both 
solitons transmit through the barrier. For Vq — 0.01, the transmitted solitons show a sign of recombining again but 
with shorter period than for the free evolution case and larger than in the case of Vq = 0.465. 

Motivated by the fact that at the coalescence point the intensity of the molecule is considerably higher than at 
other instants, we expect to find different scattering dynamics when the soliton molecule meets the interface of the 
potential at different phases of its periodic oscillation. In Fig. [TTl this is investigated by fixing the height of the 
potential barrier at Vq = 0.465 and changing the initial launching position of the molecule. Starting at Xi = 27, the 
molecule breaks after reflection. A shortly-lived molecule is obtained at Xi = 37, and a stable molecule is found for 
xi = 42, which corresponds to the Vq = 0.465 case of Fig. [TOl At xi — 52, the soliton splits at the interface into 
transmitted and reflected solitons. The coalescence point is, in this case, located at the interface. Transmission takes 
place due to the high intensity of the soliton at the coalescence point. At xi = 62, the two solitons still split as in the 
previous subfigure but with a weak transmitted soliton intensity less than 0.025 and hence will not be shown in our 
density plots which are restricted to the intensities between 0.025 and 0.15, as pointed out previously. For xi = 82, 
the coalescence point takes place before the molecule reaches the interface and both solitons reflect but the molecule 
breaks up. For xi = 92 the two reflected solitons start to recombine forming a stable molecule at xi = 102. At 
xi = 137 the reflected molecule starts to break up since the second coalescence point becomes close to the interface. 
Thus, the conclusion from this figure is that the soliton molecule is more vulnerable to break up when it meets the 
interface at the coalescence point. Equivalently, soliton molecules with larger equilibrium bond length, such that 
coalescence does not occur, will be more stable against breakup post reflection from barriers. 

Finally, we present in Figs. [T^lfni the results of scattering of a soliton molecule by a single soliton described by 
Eq. (jA3[) with normalization , center-of-mass position and speed X3 and r?3 = W3 , respectively. The effects of the 
phase, speed, and amplitude of the injected soliton are investigated separately. In Fig. [T^ . a soliton initially at 
X = —40 is launched towards a stationary soliton molecule near x = 0. At the impact, the molecule brakes up, its 
right soliton is ejected in the direction of the positive x-axis, and the left soliton combines with the scatterer soliton 
to form a new stationary molecule shifted by about a bond length to the left. We point out here that for such an 
outcome to occur, it is essential that the amplitude of the scatterer soliton is nearly equal to that of the right soliton 
of the molecule. Otherwise, a different outcome, as that of Fig. [T31 will be obtained. In Fig. [T7b . the same numerical 
experiment is repeated but with adding a tt to the phase of the scatterer soliton. Clearly, this phase addition prevents 
the formation of a new molecule resulting in three solitons diverging from each other. From applications point of 
view, the phase of tt could be used as a "key" to "unlock" the molecule for the purpose of extracting stored data. 
In Fig. [131 we show the effect of the initial speed of the injected soliton. In contrary to one's first judgment, the 
molecule preserves its structure for fast collisions, as in Fig. 113b. and breaks up for slower collisions, as in Fig. [T^. In 
Fig. [141 the injected soliton has an amplitude that is approximately two times larger than any of the two solitons of 
the molecule. The injected soliton penetrates the molecule leaving it almost unchanged apart from a center-of-mass 
shift to the left. 



IV. CONCLUSIONS 

We have used the Inverse Scattering method to derive the two solitons solution of a nonlinear Schrodiner equation 
with a parabolic potential and cubic nonlinearity with time-dependent coefficients, as given by Eq. ([2]). The solution 
was then simplified and put in a suggestive form in terms of the fundamental parameters of the two solitons, namely 
their amplitudes, center-of-mass positions and speeds, and their phases. In this form, two different regimes of the 
solution, namely the resolved solitons and overlapping solitons, were distinguished and the main features such as the 
solitons separation and relative phase were extracted. From the expression for the solitons separation we find that 
for the homogeneous case and zero solitons relative speed, the solitons separation diverges logarithmically with the 
solitons amplitude difference such that, for equal solitons, the trivial solution is obtained. 

The force of interaction was then derived, essentially, for arbitrary solitons parameters. This resulted in generalizing 
Gordon's formula 7] to i) the generalized inhomogeneous case considered here, ii) arbitrary solitons relative speed and 
amplitudes, and in) short solitons separations (compared to their width). With this formula, the possibility of forming 
soliton molecules emerged naturally, where the force at short distance was shown to be composed of an attractive 
part, resulting from the nonlinearity, and another part that oscillates between repulsive and attractive resulting from 
the time-dependent relative phase. The main features of the soliton molecule, including its equilibrium bond length 
and bond spring constant, frequency and amplitude of oscillation, effective mass, and its elastic potential energy, 
where then calculated in terms of the solitons parameters. It turns out that the amplitudes difference Ud — 112 — ni 
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plays an important role in determining these quantities. Furthermore, we show that at the condition goUs = 1, the 
two solitons coalesce while away from this condition, the solitons approach each other but remain resolved. At this 
condition, the molecule's effective mass and spring constant have maximum values. In our expressions Eqs. p2|) and 
P4p diverge because these formulae were derived assuming the solitons remain resolved. 

To have a sense of its stability we investigated numerically i) the collision of the soliton molecule with a hard wall 
and softer one, ii) scattering by a potential barrier, and in) collision with a single molecule. The first case showed 
that while the molecular structure is preserved after reflection from a hard wall, it breaks when reflecting from a 
softer one. Refiection from a finite barrier showed that the molecular structure is preserved only for specific heights 
of the barrier. For an incident molecule with the coalescence condition satisfied, the molecule will be most vulnerable 
to break up when the coalescence point takes place at the interface of the barrier. This is simply understood by 
the fact that the intensity of the soliton molecule is maximum when the two solitons coalesce such that tunneling 
becomes possible. Stability of the molecule was also investigated by scattering the molecule with a single soliton. It 
turned that slower collisions tend to break up the molecule more easily than faster ones. In addition, the outcome 
of the collision depends on the phase of the incoming soliton such that a scatterer soliton which is in-phase with the 
molecule will typically preserve its molecular structure, but for an out-of-phase soliton, the molecule breaks up. 

The two solitons solution presented here and the analysis that shows how to extract the solitons separation and 
relative phase may constitute the basis for a more accurate and detailed investigation of the origin of the soliton-soliton 
force, especially for short separations and at coalescence. The results of this paper will be hopefully of relevance to 
possible future applications of soliton molecules as data carriers or memories. 
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Appendix A: Deriving the two solitons solution 

The Lax pair associated with the Gross-Pitaevskii equation, Eq. ([2]), is obtained using our Lax pair search method 
[H ill and reads 

*^ = CJ-*-A + P-*, (Al) 
*f=iC^J-*-A-A + C(«P + a;7J)-*-A + W-*, (A2) 



where 



Mx>t) Mx>t) J' \o X2J' \o -I )' \-V9^Q* 



C{t) = %/2e''(*), Qix,t) = 5'(x,t)e(''(*)+*T'(*)'='')/2, and Ai and A2 are arbitrary constants. Here, ^'(a;,t) is the 
solution of Eq. ([2]) and $ is the auxiliary field. 

The compatibility condition — ^tx of the linear system, Eqs. (lAip and (jA2l) . is equivalent to the Gross-Pitaevskii 
equation, Eq. and its complex conjugate. For a known seed solution, \I'o(x,i), of Eq. ^ the linear system will 
have the solution $o- The Darboux transformation is defined as $[1] = $ • A — (t$, where $[1] is the transformed 
field and a = €>o • A • $o^^- Requiring the linear system to be covariant under the Darboux transformation imposes 
the transformation P[l] = P + J- (t — cr-J, where P[l] is the transformed P. This gives the new solution in terms of 
the seed solution as 

Hx,t) = M^,t) - ./A (A, _ A.)e(^(-)-^(-)--)/^ ^ . .^flf^^tf.,. (A3) 

V 50 (l)2(X,t)^pi[X,t) - (t)i{X,t)lp2[X,t) 

Using the trivial solution = as a seed, the Darboux transformation generates the well-known sech-shaped 

single soliton solution 



t) = ,/^^ e**^(-'*) sech (an(x - x„ni W)) , (A4) 
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where 

4>l{x,t) = (f)o{t) +icml(i) (a; - Xcml{t)) - ^i{t){x - Xcml{t))'^, (A5) 

Xc^i{t) = {xi e^^« + m 9{t)) e-^(*)-^°, (A6) 

g{t) = Jq e'^^'-* ''dt' , an = e^^*'^ nigQ/2, 7(0) = 70, and the constant (/)oi corresponds to an arbitrary overall phase. 
This solution corresponds to a soliton density profile that is localized at Xc■ml{i)^ moving with center-of-mass speed 
icmi(i)i and normalized to rii. 

Using this single soliton solution as a seed, the Darboux transformation generates a two solitons solution. The 
solution of the linear system, Eqs. (jAip and (|A2[) . can in this case be derived and simplified to the following form 

/exp + \ 

t) = e^+^^ ^ ^"^^^^ ^° " ^ - + 2 (2^^^l " -"""^i) > (A8) 

/ 2goni exp (-y^ + y + "^^^^^t, -2V2A| - ie-T"??! - 1 \ , , 
ib2{x,t) = e ™"i +2^^ ^ — — LL + 1 — I ii - — — , A9) 

ip,{x,t) ^ r2ix,t), (AlO) 

iP2ix,t) = -ri{x,t), (All) 

A2 = -At, (A12) 



where. 



2/ = go"ie^(*'(a;-a;cmi(t)), (A13) 
(_ie2T«g(t) (32Ai2 - g^ni^) - 16V2e-'"m^i9it) + 4.tm^g{t)) , (A14) 
(y2e-''°7n + 2iAi) Xig{t), (A15) 



y2 

ys = -AJg(i) (%/2e-T°?7i - 2iXl^ + upoi. (A16) 

Finally, the two solitons solution is obtained by substituting for ■0i.2 and 01^2 in Eq. (IA3|) . which upon substituting 
'^i — "-2 ffo/(4-\/2) + « e~'^°7'i2/ and further simplification can then be put in the form of Eq. 

Appendix B: Deriving center-of-mass positions and acceleration 
1. Center-Of-Mass Positions 

Since X and Y are functions of x, we start by examining their values near the roots of q. This task can be simplified 

by rewriting X as X = ei'="'**'5o("i(^-^<=">i(*))-"2(^-^--2(*)). poj. ^ ^ XcmiW, we get X ~ e^'''''"^3on2(x^^2it)-x,raiit)) ^ 
1, and for x ~ Xcm2{t), we have X ~ e^'^"'' 'so ni(a;cm2(t)-a:cmi(t)) Yot x > Xcm2{t) the first term in the exponent of 
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X is larger than the second one, provided that ni is not much larger than n2, as remarked at the end of this section, 
hence X ^ 1. For x < cCcmiCi) the magnitude of the first term in the exponent of X is smaller than the magnitude 
of the second one, which again leads to X ^ 1. For the region Xcmi{t) < x < Xcm2{t), the condition X ^ 1 is always 
satisfied. In conclusion, X ^ 1 for all x, apart from situations with extreme values of ni/n2. The situation is simpler 
for Y: Y ^ 1 for x > a;cm2(i), — 1 for x ~ Xcmi{t), and y ~ for x < Xcmiit)- Thus, to find Xr, we expand q in 
powers of large X and Y and to find xi and Xc, we expand q in powers of large X . 

Expanding q in powers of X and Y around oo up to first order in 1/X and 1/Y, we find 



1, , qne'^°(nd — Us)^ (qoHge'^^^'' cos z — 2ndSinz 

q = - log y - 2 log X + log - ^° ^ " ' \ — 

where 



The root of this equation gives the center-of-mass position of the right soliton. The first two log terms equal n2 go {x — 
Xcm2{t))- The third log term is constant and diverges for Ud = r/d ~ 0. The last term is small since it is proportional 
to l/X, but it is needed because it contains the phase-dependent contributions. Due to the combination of log X and 
X, finding an algebraic expression for the root of this equation will not be possible. Instead, we ignore at first the 
1 /X term to obtain the dominant contribution which will then be used to find the 1 /X contribution. This gives 

2e"'^(*) 

Xr = Xcm2{t) J ■ rlogy,-. (B3) 

Substituting back in Eq. (jBip and solving for x, we finally get 



2e ^Wlogy. 4y,".+"^ (gon,e'^"cosz,- 277rfsinz,) 
goiud + Us) go{nd + ris) {goUd'^e^^" + irjd-') 

where Zr — z{x = Xr)- 

For the left and central solitons, we expand q in powers of X; leaving Y arbitrary. In this manner, we account for 
the two roots, Xi and Xc, simultaneously. Expanding q in powers of large X, we get 

. = 2 log r + 2 log (4(y + l)^e--,.V.g + (n. + n.y)0 + ^ ' ^^^^ 

where 

C = (Y + if e-^°[2gle^^°iqdsinz{nd^ +n,:^Y)+Agoe^'> iqd^ cos z(ndY + Us) 
+ gludUsc'^^^ cos ziud + risY) + 8(y + 1)%^ sinz] 

/ {gy^ind + n^Yf+^iY+lff^d^). (B6) 
Similar to the above procedure for Xr, we solve first Eq. (jB5|) without the 1/X term, which gives 



where 



2<?2n,2e27o + Srjd^ 



(B7) 



M = gy^° {ud^ - GndUs + n^^) - 8277/. (B8) 



In the limit rid — >■ and rjd — > 0, the solution y+ approaches 0, which corresponds to xi — Xcmiit) + 
(logy+)/(ni go e''''"*'') — ^ —00. The solution y_ approaches 1, which corresponds to Xc = Xcmiit)- Thus, the so- 
lutions y^ and y_ correspond the left and central solitons, respectively. Since we will be interested only in the left 
and right solitons, we take for the rest of this section the y+. To find the contribution of the 1/X term, we substitute 
for y+ in the 1/X term of Eq. (|B5P , and then solve for Y to get a corrected expression for y+ 



Y. = 



{nd - ns)^jM+ + goe'"' e ^+ {ud ~ Ugf - 2ndn, 



-ST]d^ 

, (B9) 
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where M+ = 506^'^"'+'^^'''' (^e {ud - Us)^ - ^UdUs^ - l6T]d^ (^e + 1^ , X+ = eif«=^''("<'+"=)'=^'"+T«, Xd = 

Xcm2{t) — Xcnii{t), C'+ = C{Y — Y+,x — xi,z — zi), and zi — z{xi). Expanding for large and then solving 
for a;, we finally get 

2e-^Wlog(y+) 4^+6-^"-^^ 

go{nd-ns) X+^gle^-T" {ud^ - Qudfis + Us^) - 32r]d^ 

Final remarks about the validity of the above derivation are in order. The condition X ^ 1 will be met in the 
region x > Xcm2{t) only for ni/n2 > {Xr — Xcni2it)) / {xr — Xcmiit))- Note that the numerator of the right hand side of 
this inequality is less than the denominator by at least Xcmiit) — a;cmi(i)- For x < a;cmi(i), the condition X ^ 1 will 
be met only for ni/n2 < {xcm2{t) — xi)/{xcmiit) — xi). Here, the numerator of the right hand side of this inequality 
is larger than the denominator by at least Xcm2{t) ~ Xcmiit)- A more quantitative estimate for the ratio ni/n2 can 
be obtained using the above results for xi{ni,n2) and Xr{ni,n2)- 

Furthermore, the above-derived formula for xi is limited to values of the parameters for which the quantities M 
and M+ are positive. At M = 0, the two roots xi and Xc coincide. In Fig. [H this corresponds to the maximum of 
the g-curve lying on the x-axis. For M < 0, the maximum of the curve is below the z-axis and the two roots become 
nonreal. 



2. Acceleration 

Solitons separation A = calculated using Eqs. (jB4l) and (jBlOl) . This gives rise to Eq. (|18p with coefficients 



(n, + n,)(goWe2^°+ 477^2)' ^ 

. ^ 4.9o(y+ + 1)^ (ggnrfn,e^T°(nrf + n.y+) + Ar,d^ndY+ + n,)) 
VM{gle^'ro(^nd + nsY+)^+4{Y+ + l)^r^d^) 



^ 8{Y+ + I fr^d {g^e^^" (n/ + ns^Y+) + i{Y+ + 1)77/) 
33 = — 



Noting that: 
idW =%e^^*)-^« - .Td(t)7(t), 
id{t)=Xd{t) (7W'-7W) , 
iiit) = -ie27(*)~27o (g^^ndnse^-x' - Arjd') , 
z,(t) = - ie27(*)"2^o^(i) (52„,„,e27(°) - 4r;,2) , 

the acceleration, A(t), can be calculated to take the form of Eq. ([20]) with coefficients 



Yj,~- 
64^/1 



(B13) 



/35 = _t_^_le(-^9o^"^'=^"'-57o) (45oe2^°r,rf2 (^^/^g (^2 + 4ndn, + n^) + AP2r]d{nd + n,)) 



1 

64-v/y I 



4\ 
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64 v% 



/38 = - 
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FIG. 1: (Color online) Dashed (blue) curve: The argument of the second sech term in the two solitons solution, Eq. ((3]), 
q = a22{t){x - Xcu,2{t)) + {l/2)\og[{al + al)/{al + al)]. Thick (green) curve: The soliton intensity |*^(a;,f)|^. Light (red) 
curve: The soliton intensity with only the second sech term of Eq. ([3]). The parameters used are: 'y{t) = Q, n\ = 1, n2 = 1.01, 
x\ = 10, X2 = 20, vi = V2 = 0, go = 3, 4>oi = 4>02 = 0, and t = 0. 
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FIG. 2: (Color online) The phase of solitons (dashed curves) and solitons intensity profiles (solid curves). Long dashed curve 
(orange) corresponds to the exact phase of the solitons calculated directly from the two solitons solution, Eq. ([Sjl. Dotted 
(blue) curve corresponds to the phase of the right soliton, (fir, calculated from Eq. (|14p . Thick dashed curve corresponds to the 
phase of the central soliton, <j>c, calculated from Eq. (|15|) . Dashed-dotted curve corresponds to the phase of the left soliton, (pi, 
calculated from Eq. p6[l . The solid curves correspond to the density profiles as in Fig. [1] The parameters used are: 7(t) = 0, 
ni = 1, n2 = 1.01, xi = 10, X2 = 20, vi — 1, V2 — 1.2, go = 3, 0oi = <^02 ~ 0, and t = 0. 
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FIG. 3: (Color online) The initial acceleration of the solitons separation, A(0), versus rjd- Thick (blue) curve is calculated from 
Eq. (|20|) . Light (red) curve is calculated numerically from the exact two solitons solution, Eq. ((3]). The parameters used are: 
7(t) = 0, ni = 1, n2 — 1.01, xi — 10, X2 — 20, vi — 1, go = 3, ^oi ~ 4>02 = 0, and t — 0. 
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FIG. 4: The acceleration of the solitons separation, A(0), versus rjd and t. The parameters used are: ^{t) = 0, ni = 1, 
712 = 1.01, xi = 10, X2 = 20, vi = 1, go = 3, ^oi = 002 = 0. 
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FIG. 5: (Color online) Density profile of the soliton molecule and the solitons center-of-mass trajectories and separations, (a) 
Density plot corresponds to solitons density. Curves correspond to solitons trajectories calculated numerically from the exact 
solution Q. (b) Thick (green) curve is solitons separation calculated from the exact solution Q. Light (red) curve is the 
solitons separation calculated from formula (|25|) . (c) Density profile at some specific times. The parameters used are: 7(t) = 0, 
ni — 2.37, Ud = 0.5, xi — X2 = 10, vi = V2 = 0, go ~ 0.5, 0oi ~ 002 = 0. 





FIG. 6: Same as Fig. [5]but with m = 1.25. 
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FIG. 7: Same as Fig.[5]but with ni = 1. 
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FIG. 8: Contour plots showing the phase of the sohton molecules of Figs. 1 5171 Subfigures (a), (b), and (c), correspond to 
Fig. [5l Fig. [6l and Fig. [T] respectively. The blue (lower) and green (upper) curves correspond to the center-of-mass trajectories 
of the left and right solitons, respectively. 
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FIG. 9: Soliton molecule reflection by a potential step given by Eq. (|39|) . The dashed vertical line represents the interface of 
the potential step. The parameters used are: go = 0.5, X2 = xi = 50, (t>oi ~ 4>02 — c = 0, Ud = 0.3, ni — 1.0, rj — —0.1, 
xo = -50. (a) Vo = 100. (b),(c),(d) Vo = 0.075. 
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FIG. 10: Soliton molecule reflection by and transmitting through a potential barrier given by Eq. (I40p for different barrier 
heights. The parameters used are: go = 0.5, X2 = xi = 42, (poi = 4'02 = c = 0, Ud = 0.1, ni = 0.95, rj = —0.1, d — 0.005, 
xo = 0. 
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FIG. 11: Soliton molecule reflection by a potential barrier given by Eq. (|40p for difi'erent initial positions of the molecule. The 
parameters used are: X2 = xi, go = 0.5, (poi = </>02 = c = 0, Ud = 0.1, ni = 0.95, r] — —0.1, Vb — 0.465, d — 0.005, xo = 0. 
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FIG. 12: Collision between a single soliton and a soliton molecule for two single soliton phases, go = 0.5, X2 = xi = (poi = 
4>02 = c = 0, na = 0.2, m = 1.9, ns = 2, xs = —40, r]3 = 0.025. The phase difference between the injected soliton of (b) and 
that of (a) equals tt. 
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FIG. 13: Collision between a single soliton and a soliton molecule for two single soliton initial speeds, go = 0.5, X2 = xi = 
<poi = 002 = c = 0, nd = 0.2, m = 1, ns = 1, xs = -80. (a) 773 = -0.04, (b) 773 = -0.025. 




FIG. 14: Collision between a high intensity single soliton and a soliton molecule, go = 0.5, X2 = xi = (j>oi = ct>02 = c = 0, 
m = 0.2, ni = 1, ns = 2, 773 = 0.04, xs = -40. 



